

#include"heat_2d_eq_eq.c"


          
	m_x = (*fields_ptr).u[M_A]GRD_IND;
	m_y = (*fields_ptr).u[M_B]GRD_IND;
	m_z = (*fields_ptr).u[M_C]GRD_IND;

          psi_x = nx*auxfields.u_aux[SWAP_X_X]GRD_IND + ny*auxfields.u_aux[SWAP_X_Y]GRD_IND;
          psi_y = nx*auxfields.u_aux[SWAP_Y_X]GRD_IND + ny*auxfields.u_aux[SWAP_Y_Y]GRD_IND;
          psi_z = nx*auxfields.u_aux[SWAP_Z_X]GRD_IND + ny*auxfields.u_aux[SWAP_Z_Y]GRD_IND;

 norm_n = nx*nx*g_up_xx + ny*ny*g_up_yy + 2.*g_up_xy*nx*ny;

// if(fabs(psi_n)> 0.0001){
// printf("grid = %d, phi(%d,%d) = %f, psi(%d,%d)= %f, dx = %f, Tau = %f, delta = %f \n", grid, grid_ind1,grid_ind2,phi_n, grid_ind1,grid_ind2, psi_n, o_dx, Tau, delta);
// }
 dm_x = (m_x - m_x_n);
 dm_y = (m_y - m_y_n);
 dm_z = (m_z - m_z_n);

 dpsi_x = ((1.+ 2.*delta) + 0.0*nx + 0.0*ny)*0.5/norm_n/detg_1_2*(psi_x + psi_x_n);
 dpsi_y = ((1.+ 2.*delta) + 0.0*nx + 0.0*ny)*0.5/norm_n/detg_1_2*(psi_y + psi_y_n);
 dpsi_z = ((1.+ 2.*delta) + 0.0*nx + 0.0*ny)*0.5/norm_n/detg_1_2*(psi_z + psi_z_n);

	 (*derivs_ptr).u[M_A]GRD_IND = (*derivs_ptr).u[M_A]GRD_IND +  o_dx*(Tau*dm_x - dpsi_x); 
	 (*derivs_ptr).u[M_B]GRD_IND = (*derivs_ptr).u[M_B]GRD_IND +  o_dx*(Tau*dm_y - dpsi_y); 
	 (*derivs_ptr).u[M_C]GRD_IND = (*derivs_ptr).u[M_C]GRD_IND +  o_dx*(Tau*dm_z - dpsi_z); 

#ifdef PROJECTION

       m_x = (*fields_ptr).u[M_A]GRD_IND;
       m_y = (*fields_ptr).u[M_B]GRD_IND;
       m_z = (*fields_ptr).u[M_C]GRD_IND;

       m_t_dot_m = (*derivs_ptr).u[M_A]GRD_IND*m_x +  (*derivs_ptr).u[M_B]GRD_IND*m_y +  (*derivs_ptr).u[M_C]GRD_IND*m_z;
       norm_m = m_x*m_x + m_y*m_y + m_z*m_z;

        (*derivs_ptr).u[M_A]GRD_IND =  (*derivs_ptr).u[M_A]GRD_IND - m_x*m_t_dot_m;
        (*derivs_ptr).u[M_B]GRD_IND =  (*derivs_ptr).u[M_B]GRD_IND - m_y*m_t_dot_m;
        (*derivs_ptr).u[M_C]GRD_IND =  (*derivs_ptr).u[M_C]GRD_IND - m_z*m_t_dot_m;

#endif


// printf("grid = %d, (%d,%d) , m_x=%f, m_x_n=%f, psi_x=%f, psi_x_n=%f \n",grid, grid_ind1,grid_ind2, m_x, m_x_n, psi_x, psi_x_n); exit(0);



// Garbage


//#ifdef PROJECTION
	      //x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      //y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
                    //Normal_cub(x,y,nn_x,nn_y,nn_z,(*grid_ptr).grid);

               //m_t_nn =  (*derivs_ptr).u[M_A]GRD_IND*nn_x +  (*derivs_ptr).u[M_B]GRD_IND*nn_y +  (*derivs_ptr).u[M_C]GRD_IND*nn_z;

        //(*derivs_ptr).u[M_A]GRD_IND =  (*derivs_ptr).u[M_A]GRD_IND - nn_x*m_t_nn;
        //(*derivs_ptr).u[M_B]GRD_IND =  (*derivs_ptr).u[M_B]GRD_IND - nn_y*m_t_nn;
        //(*derivs_ptr).u[M_C]GRD_IND =  (*derivs_ptr).u[M_C]GRD_IND - nn_z*m_t_nn;

//#endif
